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1.0 SUMMARY 


A finite difference method for solving the unsteady transonic flow about harmonically 
oscillating wings is investigated. The procedure is based on separating the velocity potential into 
steady and unsteady parts and linearizing the resulting unsteady differential equation for small 
disturbances. The differential equation for the unsteady potential is linear with spatially varying 
coefficients and with the time variable eliminated by assuming harmonic motion. 

The work of this report is a direct extension of earlier studies and includes extension of a 
program using a direct solution procedure to include coordinate transformations for swept and 
tapered planforms and modification of the program to work on the new generation of vector 
computers with very large memory capabilities. 

The main results of the study are- 

1. Rederivation of the difference equations for harmonic transonic flow to include a coordinate 
transformation for swept and tapered planforms. 

2. Development of programs for three-dimensional planar configurations (including thickness) 
for the CRAY-XMP at Boeing and the CYBER VPS-32 system at the NASA Langley Research 
Center at Hampton, Virginia. 

3. An investigation of the effect of the location of the outer boundaries on accuracy for very 
small reduced frequencies. 

4. Application of the developed program to the flutter analysis of a rectangular wing of aspect 
ratio 3.0. 
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2.0 INTRODUCTION 


The development of an efficient and practical procedure for calculating unsteady transonic air 
forces for use in flutter analyses continues to be of interest to the aircraft industry. The work of this 
report extends a series of studies into a finite difference solution procedure applied to the small 
disturbance velocity potential equation, which has been linearized by assuming simple harmonic 
motion. These studies, which have been documented in References 1 through 8, have resulted in a 
pilot program that has been used successfully for the two-dimensional, typical section problem. 
Extension of the two-dimensional program to three dimensions, with the corresponding increase in 
the number of problem unknowns, resulted in a program that was very expensive to run. 

The studies described in this report include the derivation of a coordinate transformation for 
swept and tapered wings, the rewriting of a CRAY-1 program for use on the CRAY-XMP and the 
NASA Langley VPS-32 computer system in order to reduce running cost, correlation studies for 
the new pilot program for small values of reduced frequency, and application of the new program to 
a flutter example. 

The development and characteristics of the new pilot program are discussed in Section 5.0. 

The coordinate transformation, discussed in Section 6.0 and derived in the Appendix, aligns 
the mesh points with the leading and trailing edges of a planform. It is essentially a shearing 
transformation in the plane of the planform, but with additional terms in the expressions for the 
streamwise variable so that continuity of second derivatives is maintained across the planform 
edges. 

Several correlation studies are presented in Section 7.0. Generally good agreement is shown 
between pressure distributions from OPTRAN3 and a kernel function program for planforms of 
vanishing thickness. An example is presented for a wing of finite thickness. Also, it is shown that 
the boundaries of the finite difference solution region must be moved away from the wing as the 
reduced frequency is decreased in order to maintain accuracy. 

The application of the program to the flutter analysis of a wing with finite thickness is 
presented in Section 8.0. 

Finally, the authors would like to acknowledge the valuable contributions of Dr. Elizabeth L. 
Yip in the development of the solution routines, of Dr. Michael B. Bieterman in converting the 
OFTRAN3 code for the CRAY-XMP to code for the VPS-32 system, and of Dr. Robert M. Bennett to 
the calculations for the flutter example. 
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3.0 ABBREVIATIONS AND SYMBOLS 

m 

a 

Amplitude of wing oscillations 

• 

b 

Root semichord 

• 

c. 

Pressure coefficient, (p-po)/(. l / 2 p 0 U 0 2 ) where p is the local pressure, p 0 the freestream 
static pressure, and p 0 the freestream air density 


f 

Frequency in hertz 

• 

fo 

Undisturbed wing or airfoil shape 

• 

A 

Unsteady contribution to wing or airfoil shape 

A 

ij,k 

I,J,K 

x,y,Z subscripts and indices for points in the mesh 


/ 

>/=l 

• 

k 

Reduced frequency based on semichord, 2ir fb/U; same as w 


K 

Transonic parameter, (1 - M 2 )/M 2 €) 

• 

KFB 

Kernel function program for two-dimensional flow by Bland (ref. 18) 

• 

le, LE 

Leading edge 


M 

Freestream Mach number 

• 

M 0 

Mass of wing excluding support shaft 

a 

Q 

u?/e - ios(y - 1V 0joc 

^r 

A 

Qmn 

Generalized force due to pressure from the nth mode weighted with the deflection of 
the mth mode 

w 

A 

RH04 

Kernel function program for three-dimensional flow by Rowe, Redman, and Winther 
(refs. 14 and 15) 

V 

S 

Wing semispan 


t 

Time in units of bl\J 

A 

te , TE 

Trailing edge 

w 

0 
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u 

U,U 0 


x,z 

x,y,z 


Xq*Zq 

Wo* 


7 

5 

e 

h 

P 

»} 

i,n,r 

p 

PO 

p 


P<> 


Pi 

co 


w 0 


K - (7 + 1V 0jc 
Freestream velocity 

Scaled and nondimensional coordinates for two- and three-dimensional problems 
(x - X(Jb, y - yy^b, z = yZ(Jb) 

Nondimensional and unsealed coordinates for two- and three-dimensional problems 
Ratio of specific heats for air 

Jump in pressure coefficient across airfoil or wake per unit amplitude of oscillation 

Jump in at plane of wing or vortex wake 

Jump in' <fi at wing trailing edge 

Thickness ratio 

(5/M) 33 

coM/Q-M 2 ) 

Scale factor of y 0 and z 0 , y = 8^ M 33 ; also mass ratio, MoKirb^sp) 

Fraction of semispan 
Sweptwing coordinates 
Air density 

Freestream air density 

Complete, scaled perturbation velocity potential; also used for the unsteady potential 
in finite difference equations with subscripts 

Steady scaled perturbation velocity potential 

Unsteady scaled perturbation velocity potential 

Frequency in radians; also reduced frequency, same as k 

Reference frequency in radians. For the example in Figures 20 and 21, oj 0 is the 
second natural frequency of the model 
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SUBSCRIPTS 


F Flutter value 

0 Freestream value 


r 


Root or reference value 


4.0 FORMULATION AND SOLUTION 

Since the mathematical derivation of the method for the solution of the unsteady velocity 
potential for the flow about a harmonically oscillating wing is presented in Reference 1, the discus- 
sion here will be limited to a brief outline of the procedure for two dimensions. The complete 
nonlinear differential equation was simplified by assuming the flow to be a small perturbation 
from a uniform stream near the speed of sound. The resulting equation for unsteady flow is 

[K — (7 — 1) <P, - (7 + 1) - (Pxx + <Pyy - (2 <p xt + *?„)/ e = 0 (!) 

where K = (1-M 2 ) / (M 2 «), M is the freestream Mach number of velocity U 0 in the ^-direction, 
x and y are made dimensionless to the semichord b of the airfoil and the time t to the ratio 
b/ U 0 . With the airfoil shape as a function of time defined by the relation 

y 0 = 5f{x,t) 


the linearized boundary condition becomes 


<fiy=f x (X,t)+f,(X,t) 


(2) 


The quantity 5 is associated with properties of the airfoil (such as maximum thickness ratio, 
camber, or maximum angle of attack) and is assumed to be small. The coordinate y is scaled to the 
dimensionless physical coordinate y 0 according to 

y = 6 1/3 M 273 ;^ 


and e is given in terms of 5 by 


« = (5/M) 273 

The pressure coefficient is found from the relation 

Cp = -2e (<?,+ <?,) 


The preceding differential equation is simplified by assuming the velocity potential to be 
separable into a steady-state potential and a potential representing the unsteady effects. We write 
for the perturbation velocity potential 

v = -<p 0 (x,y) + <Pi (x,y)e* u ‘ (3) 


and for the body shape 

y 0 = bf(x,t) = 5 [/„(*) + f x (*)*'“']. 
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Since the steady-state terms must satisfy the boundary conditions and the differential equa- 
tion in the absence of oscillations, we obtain 

[K - (7 + 1 ) <p 0x ] <p 0xx + <p 0yy = 0 (4) 

with 

<Po y = fo(x), y = 0, -1<*<1 (5) 


On the assumption that the oscillations are small and products of may be neglected, equa- 
tions (4.1) and (4.2) with the aid of equations (4.4) and (4.5) yield 

[[K — (7 + + — (2/w/e) <P\ X + Q<P\ = 0 (6) 

where 

q = w 2 /e - /oj (7 - l)<p 0xx 
subject to the wing boundary conditions 

<Pi y = fi x + 'Vi (x), y = 0, -1<*<1 (7) 

The boundary condition that the pressure be continuous across the wake from the trailing 
edge was found in terms of the jump in potential A <fi x to be 

M = A Vli e- iu (x - x “> (8) ' 

where is the jump in the potential at x=x, e just downstream of the trailing edge and is 
determined to satisfy the Kutta condition that the jump in pressure vanish at the trailing edge. 
The quantity A^ is also used in the difference formulation for the derivative ^ , ) v to satisfy conti- 
nuity of normal flow across the trailing edge wake. 

For the set of difference equations to be determinate, the boundary conditions on the outer 
edges of the mesh must be specified. In the original unsteady formulation, these boundary condi- 
tions were derived from asymptotic integral relations in a manner parallel to that used by Klunker 
(ref. 9) for steady flow. A later formulation (ref. 3) applies an outgoing plane wave boundary 
condition to the outer edges of the mesh. This boundary condition is numerically simpler to apply 
and is equivalent to the first-order nonreflecting boundary conditions derived by Engquist and 
Majda (ref. 10). 

A computer program for solving the steady-state transonic flow about lifting airfoils based on 
equations (4.4) and (4.6) was developed by Cole, Murman, and Krupp (refs. 11 
and 12). 
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The similarity of the unsteady differential equation to the steady-state equation suggests that 
the method of Cole, Murman, and Krupp should be an effective way to solve equation (4.6) for the 
unsteady potential Note that equation (4.6) is of mixed type, being elliptic or hyperbolic when- 
ever equation (4.4) is elliptic or hyperbolic. Central differencing was used at all points for the y 
derivative and all subsonic or elliptic points for the x derivatives. Backward (or upstream) differ- 
ences were used for the x derivatives at all hyperbolic points. The preferred numerical approach to 
solving the resulting large-order set of difference equations is a relaxation procedure, which per- 
mits the calculation to be made as a sequence of relatively small problems. However, as discussed 
in preceding NASA reports by the authors (refs. 3 and 4), a significant problem of solution conver- 
gence with the relaxation procedure was encountered that severely limits the range of Mach num- 
ber and reduced frequency for which solutions may be obtained. Accordingly, an out-of-core solver 
(ref. 8) was developed to solve the complete set of difference equations simultaneously, which for 
two-dimensional flow is relatively efficient. 

The size of practical three-dimensional problems is such that the program of Reference 8 was 
too expensive for practical application. In Reference 8, the authors explored alternatives to the 
direct solution procedure and ways of making the direct solution more efficient. The conclusions 
from this work provided the motivation for work discussed in the following 
sections. 
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5.0 THE PILOT PROGRAM FOR LIFTING SURFACES 


The initial block of work was to develop a pilot program for three-dimensional lifting surfaces 
utilizing the harmonic procedures for unsteady transonic flow developed in References 1 through 8. 
This program was to be designed for use on the CRAY-XMP and was to be a refinement of a pilot 
program developed for the the CRAY-1. This initial program was developed and tested for a 
machine with 1 million words of core storage and was written during the installation at Boeing of 
both the CRAY hardware and software. The Boeing CRAY-XMP has 4 million words of core storage, 
which permits working with larger blocks of data in memory. Larger blocks mean a smaller 
number of longer do-loops; vectorizing the longer do-loops reduces execution time. Another key 
element of the XMP is its solid-state disk (SSD) with some 128 million words of storage and a data 
transfer rate 100 times that of a conventional disk. These were significant elements in the 
motivation to refine a program for which I/O requirements accounted for 75% of the run cost. 

After completion, the CRAY-XMP version of the pilot program was modified for use on the 
VPS-32 system at NASA Langley. The CRAY-XMP version, for example, included calls to a number 
of subroutines in VectorPak and BCSLIB. VectorPak is a subroutine library of carefully optimized 
mathematical subprograms for the CRAY-XMP and some of the other vector computers. BCSLIB is 
a comprehensive, general-purpose library of mathematical, statistical, and utility subprograms. 
Since BCSLIB and VectorPak are not available on the NASA Langley VPS-32 system, OPTRAN3 
was augmented with FORTRAN equivalents of the required subroutines. Conversion was 
accomplished in a relatively direct fashion, and the resulting program is not as efficient on the 
VPS-32 system as it is on the CRAY-XMP. Relative run times are discussed below in Section 5.2. 

5.1 PROGRAM CHARACTERISTICS 

The program is set up for a single, noninterfering surface with straight leading and trailing 
edges. It includes a coordinate transformation to align the computational mesh with the leading 
and trailing edges of swept and tapered planforms. It includes capabilities for calculating the 
downwash for vertical translation, pitch, flapping, or full span control surface modes. However, 
with relatively simple changes in the subroutine that calculates the downwash distribution 
(subroutine AIRFOIL), the program can be adapted to run user-supplied modal deflections and 
slopes (as was done for the flutter calculation in sec. 8.0). The modal information is supplied to the 
program only at the mesh point locations. No specific information is provided for discontinuities in 
slopes in the deflection pattern such as occur in (for example) a control surface mode. However, if 
mesh points are concentrated about the hingeline, the proper gradients in the pressure pattern will 
be observed. The program has no means for calculating slopes from input deflections. 

Program output includes pressure distributions and sectional and generalized forces. 
Relatively simple changes to the subroutine CPR would permit obtaining the pressure or force 
data in nearly any form desired. 

Care must be taken in calculating sectional generalized forces from pressure distributions 
defined with a limited number of finite difference points. For configurations with subsonic leading 
edges and thus pressure distributions that are singular at the leading edge, use of trapezoidal 
integration may significantly underestimate the total sectional force. In OPTRAN3, an assumed 
pressure distribution with a square root singularity is used in the vicinity of the leading edge. The 
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second mesh point aft of the leading edge, x(I0 + 1), is used as a critical point. The force due to the 
pressure ahead of this point is obtained analytically using the formula 


AC „ = AC, (*(I0 + 1)) 


» Vx(I0 + 1) + 1 
Vx + 1 


The force due to the pressure aft of x(I0 + 1) is obtained using trapezoidal integration. Using 
this procedure, results from OPTRAN3 correlated well with results from RH04 (see sec. 8). 
Pressure and sectional force data were not available for correlation with a finite thickness case. 
Care should be taken in applying this procedure to thick wings with blunt leading edges for which 
the usual rapid rise in pressure at the leading edge may be attenuated. 

For purposes of this report, the number of points in a mesh is given in terms of a full mesh 
about the right-hand side of the planform. Thus, an IMAX x JMAX x KMAX mesh has IMAX 
points in the x-direction (streamwise) with x(l) on the upstream boundary and x(IMAX) on the 
downstream boundary, JMAX points in the y-direction (spanwise) with y(2) on the planform 
centerline, y( 3) = -yd), and y(JMAX) on the outboard boundary, and KMAX points in the 
z-direction (vertical) with z(l) on the lower boundary and z(KMAX) on the upper boundary. A 
full-mesh solution for problems without symmetry about the x-y plane has 
(IMAX - 2) x (JMAX - 2) x (KMAX - 2) unknowns. A half-mesh solution, for problems with 
symmetry about the x-y plane and for which only points below the planform need be included, has 
(IMAX - 2) x (JMAX - 2) x (KMAX/2 - 1) unknowns. 

In similar fashion, for two-dimensional problems, an IMAX x KMAX mesh would have 
(IMAX - 2) x (KMAX -2) unknowns for a full solution and (IMAX - 2) x (KMAX/2 -1) unknowns 
for a half-mesh or symmetric solution. 

All the examples in this report are problems with symmetry about the x-y plane and thus are 
half-mesh solutions. 

5.2 PROGRAM RUNNING TIME 


The new program, when run on the CRAY-XMP, is significantly more efficient than the CRAY-1 
version. The cost of a typical run for the new program in CPU seconds is one-tenth that of the 
original program. Also, due to the use of the SSD, the I/O cost for the refined program is 25% of the 
total run cost according to the Boeing cost algorithm rather than the 75% for the CRAY-1 program. 

The running time for the new program on the CRAY-XMP for a half-mesh solution using 
48 x 17 x 32 mesh (10,350 unknowns) is about 120 CPU seconds. The number of right-hand sides 
(the number of mode shapes) has a small effect on the running time. The running time increases 
roughly as a cubic of the increase in the number of spanwise and/or vertical points. The running 
time variation is roughly linear with increasing points in the streamwise direction. A plot of 
running time versus the number of points spanwise is shown in Figure 1, which generally confirms 
these variations in running times. The cubic curve is based on JMAX -2, where JMAX is defined 
as the number of spanwise points in the mesh including one point on the wing centerline and one 
point to the left of the centerline. 

The VPS-32 version of the program requires about ten times more CPU seconds than the 
version on the CRAY-XMP. 
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6.0 A TRANSFORMATION FOR SWEPT AND TAPERED PLANFORMS 


A very desirable characteristic of finite difference programs for three-dimensional lifting 
surfaces is the matching of the calculation mesh to the planform geometry of swept and tapered 
wings. This is particularly important in the vicinity of surface leading edges, where there are large 
gradients in the pressure distributions, and for highly tapered wings, where, with rectangular 
grids, it is difficult to obtain a practical number of points along both the root and tip chords. 

The coordinate transformation derived in the Appendix and included in the pilot program 
generally aligns mesh points with the leading and trailing edges of a swept and tapered planform 
using a stretching and shearing algorithm. The calculation volume is divided into three regions by 
vertical planes through the leading and trailing edges of the planform. These extend from the 
centerline to the side boundary. The upstream and downstream boundaries are perpendicular to 
the root plane. 

The transformation is set up for wings with straight leading and trailing edges. For wings with 
swept edges, the apex at the centerline is rounded to remove singularities in the distributions. The 
rounding is accomplished so that the new edges have continuous second-order derivatives with 
respect to the spanwise variable. 

The swept edges are also extended to meet the outside boundary at right angles. Here again 
the functions are made continuous in the second derivative with respect to the spanwise variable. 

Generally, this transformation conforms to the suggestion of Goorjian and Guruswamy in 
Reference 13 and used in the XTRAN3 programs. The transformations used here for the upstream 
and downstream regions employ additional nonlinear terms in the expressions for the streamwise 
variable. 

A complete derivation of the transformation and its incorporation into the finite difference 
equations is presented in the Appendix. Section A.l presents the basic transformation equations. 
Implementation requires definition of the planform leading and trailing edges from the wingtip to 
the outer (side) mesh boundary. This may be accomplished in an automated fashion using the 
equations presented in A.2, which are included in OPTRAN3. However, this procedure will work 
well only if the planform has a moderate taper ratio. The problem is that the location of the outer 
boundary is determined by intersection of extensions of the edges. Thus if the taper ratio is large, 
the outer boundary will be too close to the wingtip; conversely, if the taper ratio is small, the outer 
boundary will be too far from the wingtip. Also, the procedure of A.2 assumes the planform edges 
are straight and without breaks. As an alternative, the leading and trailing edge extensions may 
be worked out by hand and entered as data input. For OPTRAN3, this means modifying 
subroutine PLNFORM to read the edge extension data rather than calculating it. The 
transformation of A.l is applied to the equation for harmonic transonic flow in A. 3.1, and the 
resulting finite difference equations are derived in A.3.2 and A.3.3. 
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7.0 CORRELATION OF PRESSURE DISTRIBUTIONS 


Validation of the new pilot program was accomplished by matching results from the new pilot 
program for OPTRAN3 with distributions from a reliable kernel function program, RH04 (refs. 14 
and 15). Discussions in this section will center mainly on correlation with RH04. For planforms of 
vanishing thickness and thus linear flow, results from OPTRAN3 and RH04 should match 
exactly. 

Section 7.1 presents an example for a swept and tapered planform. Section 7.2 discusses the 
effects of the location of the outer boundaries for small values of reduced frequency. 

7.1 A SWEPT AND TAPERED WING EXAMPLE 

Pressure distributions are presented in this section for a swept and tapered wing of vanishing 
thickness. A sketch of the planform is shown in Figure 2. The surface studied has a 30-deg swept 
leading edge and a taper ratio of 0.5, with a semispan equal to 1.25 semichords. The calculations 
were made for a Mach number of 0.8 and with the wing oscillating in pitch about the root midchord 
at a reduced frequency of 0.06 based on the root semichord. The results should match pressure 
distributions from a linear solution. For this study, results are correlated with distributions from 
the RH04 program described in references 14 and 15. 

Results from OPTRAN3 correlate well with corresponding results from RH04. Examples of 
this correlation for the root and tip (j; = 0.8) chords are shown in Figures 3a through 4b. 
Distributions are shown from RH04, and from OPTRAN3 for half-mesh solutions using a 
30 x 12 x 20 mesh and a 44 x 14 x 32 mesh. For the root chord, all three calculations include 
rounding of the leading and trailing edge. For the real part, there is little difference between the 
two OPTRAN3 calculations and RH04. Both meshes do a good job of predicting the leading edge 
singularity. For the imaginary part, the finer mesh does a noticeably better job of matching the 
RH04 results over the the forward part of the wing. For the 80% spanwise chord, the finer mesh 
again does a better job of matching RH04. Although the difference is marginal for the real part, it 
is significant over the leading edge region for the imaginary part. 

For both meshes, the upstream and downstream boundaries were at x = - 2.0 and 3.0 and the 
side boundary was at y = 4.5. The upper and lower boundaries were at z = + 4.5 for the fine mesh 
and z = +3.0 for the coarse mesh.* In view of the subsequent difficulties with the location of the 
mesh boundaries at low reduced frequencies, results might well be improved by moving the 
boundaries further out. This point will be discussed in the next section. 

7.2 A FINITE THICKNESS EXAMPLE 

Pressure distributions are presented in this section for a rectangular wing with a 10% thick 
parabolic arc airfoil section. The semispan is 1.5 times the semichord, and results are presented for 


* The x, y, and z values are given here in units of root semichord. Thus the upstream boundary 
is 1/2 root chord length in front of the upstream boundary at the root; the downstream boundary is 
one root chord length aft of the trailing edge; and the spanwise boundary is four semispans from 
the root. 
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M = 0.9 and k = 0.06. The mesh used is 44 x 12 x 32 with the boundaries relatively close to the 
wing: -2.0 s r s 3.0, y < 4.5, -4.5 < z < 4.5. The flow condition is symmetric (a = 0) and the 
wing is oscillating in pitch about the leading edge. The steady velocity potential distribution was 
calculated with XTRAN3S (ref. 16). Computer resources required are essentially the same as 
required for the flat plate case. Results are presented in Figures 4a and 4b for the root section and 
for the 80% spanwise station. We do not have matching results from another source, and thus the 
most that we can say is that the results appear reasonable. We have run the OPTRAN3-ADI 
program (OPTRAN3 with an alternating direction solution package, ref. 8) with essentially the 
same results. It should be noted that other three-dimensional correlation studies are presented in 
References 4 through 8. It is surprising to find a blip in the curves of Figures 4a and 4b at the sonic 
line, where the flow is changing from subsonic to supersonic flow on the forward half of the section. 
Changing the Mach number and the point distributions (slightly in both cases) made no difference 
in this blip. On looking back at the calculations in Reference 8, hints of a blip appear in the 
two-dimensional results for the NACA-64A010 airfoil (see figs. 8 and 9 of ref. 8) but not in the 
three-dimensional results for the circular arc airfoil. 

7.3 SENSITIVITIES TO BOUNDARY LOCATIONS 

Initial results for the flutter analysis indicated a problem with the generalized forces for small 
values of reduced frequency, in this particular case for a reduced frequency of 0.01 based on the 
root semichord. It is possible that this is really a problem of Mach number as well as reduced 
frequency since the solution procedure was found to be sensitive to the parameter 

X = wM/(l - M 2 ) 

In this study, only sensitivities to reduced frequency have been examined. The concern was that 
the imaginary parts of the generalized forces were not going to zero as the reduced frequency went 
to zero. Indeed, correlation of the imaginary part of the pressure between OPTRAN3 and RH04 
was very poor at low values of k. 

The problem appeared during the flutter analysis of the wing in Section 8.0 and may most 
easily be illustrated by plotting the imaginary part of the Q 2 2 term (the generalized force due to the 
pressure from the second mode weighted with itself) versus reduced frequency. The second elastic 
mode of this wing, a rectangular wing with an aspect ratio of 3.0, is essentially a pitch mode. The 
imaginary part of the Q 22 term should go to zero as k goes to zero. Results for the rectangular wing 
of aspect ratio 3.0 and of vanishing thickness are presented in Figure 5 for both OPTRAN3 and a 
reliable three-dimensional kernel function program, RH04. It is seen that the results from RH04 
do tend to zero as k goes to zero, while those from OPTRAN3 as calculated with a relatively small 
finite difference volume do not. However, even at k = 0.1 the two imaginary parts of Q 22 are 
relatively close together and, although not shown in Figure 5, become even closer at higher values 
of reduced frequency. 

For reference purposes, the pressure distributions for the second mode for M = 0.8 and k = 0.3 
from OPTRAN3 and RH04 are shown in Figures 6a through 6d. OPTRAN3 results are shown for 
two different meshes, both with a 44 x 14 x 32 mesh and both with the side boundary at y = 4.5: a 
“small” mesh with the upstream boundary at x = - 2.0, the downstream boundary at x = 3.0, and 
the upper/lower boundaries at z - +4.5, and a “large” mesh with the upstream/downstream 
boundaries at jc = + 8.0 and the upperAower boundaries at z = + 8.0. Results are presented for two 
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chords, the root chord and a chord near the tip. First, for this particular case, the correlation 
between OPTRAN3 and RH04 is good for both the “small” and the “large” meshes and for both 
chords. Moreover, the two OPTRAN3 results are noticeably more like each other than either of 
them is like the RH04 result. Thus, in this case, the larger solution area does not affect the result. 

The phenomenon of deteriorating accuracy as the reduced frequency went to zero had been 
noted in the two-dimensional problem by Siedel, Bennett, and Whitlow (ref. 17) and was resolved 
by moving the upstream boundary further upstream. This was immediately tried on the 
three-dimensional problem but with only limited success. It ultimately proved that all the mesh 
boundaries had to be moved out to obtain the desired correlation between OPTRAN3 and RH04. 

The next part of the discussion centers on results for a two-dimensional airfoil section of 
vanishing thickness oscillating in pitch about the leading edge. The results will be presented in 
terms of the pressure distributions. It is recognized that, for small k values, the amplitudes of the 
imaginary parts of the generalized force are small relative to the real part. Indeed, the difference 
in imaginary parts would not be noticeable if they were plotted on the same scale as the real parts. 
However, in light of the importance of the imaginary part of Q 22 having the proper behavior as k 
goes to zero, it is of interest to investigate this problem in some detail. Due to the computer costs 
entailed in running the three-dimensional program, it was decided to begin by comparing results 
from OPTRAN2 and a well-tested two-dimensional kernel function program by Bland (ref. 18), 
hereafter called KFB. 

Figures 7 through 10 present results for M = 0.8 and k values of 0.3 to 0.001 for a 
“reasonable-sized” two-dimensional mesh: a 72x60 point pattern (see sec. 5.1 for definition of 
mesh sizes) with |x| < 10 and jy| < 10. The intent is to show (1) the relative amplitudes of the real 
and imaginary parts and how they change with k value, and (2) how the correlation between 
OPTRAN2 and KFB degrades as the Ar-value is reduced while the solution area is held constant. 
The correlation at k = 0.3 (fig. 7) is excellent. The falloff in correlation at k = 0.1 (fig. 8) is barely 
discernible. At k - 0.01 (fig. 9), the real parts still coincide, but there is a noticeable difference in 
the imaginary parts. This is the expected behavior where, for a given finite difference mesh sized 
to provide good answers at nominal values of k, accuracy is lost at very small values of k. However, 
at k = 0.001 (fig. 10), the pressures from the two methods appear to again coincide. This latter 
result is considered fortuitous for reasons illustrated in Figure 11, which presents plots of the 
imaginary part of the pressure for several upper/lower boundary locations. In this case, the 
correlation between OPTRAN2 and the KFB method decreases as the boundary is moved out, 
although not in the expected sequence. Here, the correlation is best with the upper/lower 
boundaries at y = + 10, significantly worse with the boundaries a y - +50, and then somewhat 
better with the boundaries at y - + 100. This inconsistency is probably due to the need for more 
mesh points as the boundaries are moved out. Finally, it is noted that in this particular case 
(k - 0.001) the imaginary parts are very small and the differences discussed may not be very 
significant. 

An example of the pressure distribution from OPTRAN2 moving across the results from the 
KFB method is shown in Figure 12 for k = 0.01. Here, the boundaries are moved from y = + 10 
(the OPTRAN2 result is below that from the KFB method) to y = +30 (the OPTRAN2 result is 
above that from the KFB method), to y = + 100 (little change in the OPTRAN2 result). When the 
upstream/downstream boundaries are set at x - +40 and with the upper/lower boundaries at 
y - + 60, the OPTRAN2 results correlate with the KFB results very well. 
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A three-dimensional example is presented in Figures 13a through 14b for a full-span, aspect 
ratio 3.0 rectangular wing of vanishing thickness and oscillating in pitch about the leading edge. 
Pressure distributions for OPTRAN3 for four different mesh areas and for RH04 are presented for 
two chords, the root chord and for rj = 0.885. Figures 13a and 13b present both the real and 
imaginary parts of the pressure, showing the relative amplitudes of the two parts and with very 
good correlation between all five sets of results for the real part. The RH04 results, represented by 
the solid line, are “buried” in the curves for the three large-area solutions for the real part, and are 
the uppermost curve for the imaginary part. This latter curve may be seen more clearly in Figures 
14a and 14b, which show the imaginary part of the pressure to a larger scale. As the outboard 
boundary is moved out by adding mesh points (the original spanwise point distribution is left 
unchanged) the imaginary part of the pressure moves toward the RH04 results. 

The reasons for needing to move the boundaries out for very small values of reduced frequency 
is not completely clear. It would seem logical that the boundaries must be far enough out and 
enough mesh points must be used so that the wave pattern in the potential distribution would be 
well defined. Samples of the two-dimensional velocity potential for k = 0.001 along a constant 
x-line are shown in Figures 15a and 15b for a column of points just aft of the leading edge. The 
pressure is calculated by extrapolating the velocity potential to the wing surface ( z = 0) along each 
column of points and then calculating the pressure coefficient as 

C p = 2 e (<p x + <p,) 

The pattern is shown out to z = 16 in Figure 15a and then an enlarged section - 2.0 < z < 0 is 
shown for Figure 15b. Although the patterns do vary as the upper/lower boundaries are moved 
away from the airfoil, the differences over the first two points (which are the points used in the 
interpolation) are very subtle. One would not expect to see the differences in the pressure 
distributions of Figure 11, which are calculated from the potential distributions of Figure 15b. 

Generally, these examples show that, for small values of reduced frequency, the outer 
boundaries of the solution mesh must be much further from the wing planform than for larger 
values of k. There will always be a dilemma of using enough mesh points to properly define the 
waves in the velocity potential field while keeping the number of points to a minimum to reduce 
computer costs. It is practical to vary the location of the outer boundaries with k value. However, it 
is recommended that the meshes be evaluated by correlating results from OPTRAN3 and a kernel 
function program such as RH04. 

In conclusion, compared to higher values, lower values of reduced frequency require more 
remote boundaries and larger numbers of mesh points, resulting in larger computer costs, in order 
to retain accuracy. 



8.0 A FLUTTER EXAMPLE 


OPTRAN3 was applied to the flutter analysis of the rectangular wing studied by Cole in 
Reference 19. This wing, with an aspect ratio of 1.5 and a 64A010 airfoil section, had been tested in 
the LaRC 0.3m Cryogenic Wind Tunnel (TCT). The model consisted of a rigid wing with an integral 
flexible support that was cantilevered from the wind tunnel wall. Experimental flutter points were 
obtained over the Mach number range of 0.5 to 0.9. There were no signs of transonic flow effects on 
the flutter results at the higher Mach numbers. Analytical results were obtained using the 
NASA-Langley Flutter-Analysis-System (ref. 20), FAST, which uses subsonic lifting surface theory 
analogous to RH04 (ref. 21) to calculate the aerodynamic forces. The flutter analysis was 
performed using a combination of calculated and measured properties. A set of four natural modes 
was used as generalized coordinates. The FAST results “. . . predicted the experimental flutter 
dynamic pressures rather well, with analytical results ranging up to 5 percent nonconservative at 
the lowest Mach number tested and 5 percent conservative at the higher Mach numbers tested” 
(ref. 19). 

To demonstrate its application of OPTRAN3 to the flutter problem, OPTRAN3 was used to 
calculate a flutter point at M = 0.8. The calculations were made on the NASA-Langley VPS-32 
computer system using NASA-supplied modal and velocity potential data as input. The resulting 
generalized forces were used in FAST for a V-g type flutter solution. 

The flutter analyses were made using the first four natural modes as generalized coordinates. 
The modal deflections and slopes were interpolated from the NASA structural grid to the 
aerodynamic control points used by OPTRAN3. A steady state velocity potential distribution, 
calculated with the NASA-Langley CAP-TSD program (refs. 22 and 23), was also interpolated for 
input. The effects of a thickness distribution are included in OPTRAN3 through this steady 
potential, which is used in calculating the spatially varying coefficients of the finite difference 
equations. For the current problem, a set of generalized forces for five A:-values over the range of 
0.01 to 0.2 with forces for intermediate ^-values being found by interpolation. A separate 
OPTRAN3 run is required for each Ar-value. However, each run includes all the mode shapes, and 
additional modes do not significantly effect the cost of a run. Computer costs, then, are linear with 
respect to the number of reduced frequencies for which generalized forces must be calculated, but 
are nearly independent of the number of generalized coordinates used. 

The steady pressure distribution calculated from the velocity potential used as input to 
OPTRAN3 is presented in Figures 16 and 17. The steady chordwise pressure distribution at the 
wing root shown in Figure 16 has a small supersonic region centered about a point 40% of the 
chord aft from the leading edge, and there is no discernible shock at the aft end of the supersonic 
region. Thus the flow is nearly subsonic at the root, and this pattern holds outboard as shown by 
the distribution for the complete planform in Figure 17. 

Figures 18a through 18d present the generalized air forces for the first two modes. Here, the 
kernel function forces from FAST are compared with forces from OPTRAN3 both with and without 
thickness. Generally, the three curves lie very close together for all four generalized forces and for 
both the real and imaginary parts. Differences between the two flat plate cases (FAST and 
OPTRAN3 without thickness), and small inconsistencies between the flat plate and with-thickness 
results, such as Q 12 (fig. 17) where FAST and OPTRAN3 with-thickness results match exactly and 
the forces from OPTRAN3 without thickness are slightly different, are probably caused by the 
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limited number of mesh points used as a result of practical restrictions on available computer 
resources. In this case a relatively coarse grid of 48x17x32 was used. In addition, due to the low 
/r-values at which flutter was expected, it was important to have the boundaries of the finite 
difference mesh as far from the planform as possible. After consideration of the number of mesh 
points to be used, the grid boundaries were set at 

-10.0 < x < +10.0 y < 12.75 -10.0 :£ z < +10.0 

Improved correlation with the FAST results would be expected from using more mesh points for the 
existing solution volume. Even better correlation would be expected if, with the increased number 
of mesh points, the boundaries were moved out even further. 

As would be expected from looking at the generalized forces, the flutter results for the three 
calculations also match each other very closely. A large variation in the flutter-reduced 
frequency-with-mass ratio (defined as M 0 /(it6 2 j p)) is shown in Figure 19, with the critical lvalue 
for the problem being solved here of order 0.05. Figures 20 and 21 present the flutter speed index 
(defined as U F /(h w and flutter frequency ratio (defined as u F /w) versus mass ratio. Both the 
flutter speed index and the frequency ratio are relatively insensitive to the mass ratio. The 
exception is the flutter speed index for low values of reduced frequency, which shows the expected 
rapid increase in index at very small values of mass ratio. 

The curves of flutter speed index versus mass ratio are smooth and essentially parallel to each 
other. Assuming the FAST results match the experimental results, OPTRAN3 without thickness 
provides a flutter speed index curve that is between 2% and 3% nonconservative, and OPTRAN3 
with thickness provides a boundary that is between 1% and 2% nonconservative. 
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9.0 CONCLUSIONS 


The work for this report includes a derivation of a coordinate transformation for swept and 
tapered wings, reprogramming a CRAY-1 program for the CRAYXMP and the Langley VPS-32 
computer system, an investigation of the effect on accuracy of the location of the outer mesh 
boundaries, and a sample flutter calculation of a low-aspect-ratio rectangular wing. 

The coordinate transformation for swept and tapered wings is essentially a shearing 
transformation in the plane of the planform, but with additional terms in the expressions for the 
streamwise variable so that continuity of second derivatives is maintained across the planform 
edges. This transformation was included in a pilot program and checked out on a planform of 
vanishing thickness. 

An existing CRAY-1 program was rewritten to include the transformation described above for 
the Boeing CRAY-XMP and the NASA Langley VPS-32 computer system. The CRAY-XMP version, 
which made full use of the solid-state disk and the larger available core storage, proved to be 
nearly ten times more efficient than the CRAY-1 program. The VPS-32 version of the program was 
not nearly as efficient as the CRAY-XMP program. 

Even with these improvements in efficiencies, the CRAY-XMP program is marginally efficient 
for practical flutter evaluation problems. However, another version of the program developed by 
Rowe and Ehlers (ref. 21), in which improved boundary conditions permit a reduction in mesh size 
and (ultimately) in the number of unknowns, may provide a means for reducing computer costs to 
the point where these harmonic finite difference programs would be practical. 

A study was made of the effect of mesh boundary location on accuracy for small values of 
reduced frequency. It was found that for smaller A:-values, the boundaries had to be moved further 
out to retain accuracy. Relatively good correlation was obtained down to k = 0.01 at M = 0.8. 

Finally, the application of OPTRAN3 to a flutter analysis was demonstrated with calculations 
for a rectangular wing. 
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Figure 2. Planform for Swept and Tapered Wing Example 
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Figure 3a. Comparison of Pressure Distributions From 0PTRAN3 and RH04 for a Swept 
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Figure 3b. Comparison of Pressure Distributions From OPTRAN3 and RH04 for a Swept 
and Tapered Wing of Vanishing Thickness; M = 0.8, k = 0.06, Root Chord, 
Imaginary Part 


29 






RH04 


AC p , 

imaginary 


Figure 3d. 


M = 0.8, k = 0.06 
i = 0.80 

Pitch about root midchord OPTRAN3; 

-2.0 ^ x ^ 3.0 
0 ^ y £ 4.5 
|Z| ^ 4.5 


0 44x14x32 

X 30x20x12 



Comparison of Pressure Distributions From OPTRAN3 and RH04 for a Swept 
and Tapered Wing of Vanishing Thickness; M = 0.8, k = 0.06, rj = 0.8, Imaginary 
Part 


31 


M = 0.9, k = 0.06 
Aspect ratio = 1 .25 
Pitch about leading edge 
5 = 0.10, parabolic arc airfoil section 
44x14x32 mesh 
-2.0 ^ x ^ 3.0 
0 <; y ^ 4.5 
-4.5 ^ z ^ 4.5 


+ Root chord 
0 ij - 0.8 


-2i — 
- 1.0 





Figure 4a. Pressure Distributions for a Rectangular Wing With a Parabolic Arc Airfoil 
Section; M = 0.9, k = 0.06, Real Part 


+ Root chord 
0 n = 0.8 


M = 0.9, k = 0.06 
Aspect ratio = 1 .25 
Pitch about leading edge 
5 = 0.10, parabolic arc airfoil section 
44x14x32 mesh 
-2.0 ^ x ^ 3.0 
0 £ y £ 4.5 
-4.5 ^ z ^ 4.5 



x 

Figure 4b. Pressure Distributions for a Rectangular Wing With a Parabolic Arc Airfoil 
Section; M = 0.9, k = 0.06, Imaginary Part 



k-Value 


Q 22 , 

imaginary 



Figure 5. Imaginary Part of the Total Generalized Force Versus k-Value for NASA Mode 2 for 
M = 0.8 
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Rectangular wing, aspect ratio = 3.0 
Pitch about leading edge 

RH04 

OPTRAN3; 44x14x32, y £ 4.5 
□ |x| £ 8.0, |z| ^ 8.0 

* -2.0 ^ x ^ 3.0, |z| <, 4.5 


-10 4 - 
- 1.0 



Figure 6b. Comparison of Pressure Distributions From OPTRAN3 for Two Outer Boundary 
Locations and for RH04; Rectangular Wing, M = 0.8, k = 0.3, Root Chord, 
Imaginary Part 



Rectangular wing, aspect ratio = 3.0 
Pitch about leading edge 

— RH04 


OPTRAN3; 44x14x32, y ^ 4.5 
□ |x| ^ 8.0, |z| ^ 8.0 

X -2.0 ^ x ^ 3.0, |z| ^ 4.5 
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Figure 6c. Comparison of Pressure Distributions From OPTRAN3 for Two Outer Boundary 
Locations and for RH04; Rectangular Wing, M = 0.8, k = 0.3, ij = 0.885, Real 

Part 
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Figure 6d. Comparison of Pressure Distributions From 0PTRAN3 for Two Outer Boundary 
Locations and for RH04; Rectangular Wing, M = 0.8, k = 0.3, rj = 0.885, 
Imaginary Part 
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M = 0.8, k = 0.01 
Pitch about leading edge 
+ KFB 

0 OPTRAN2; 72x60, |x| £ 10, |z| £ 10 



x 


Figure 9. Comparison of Pressure Distributions From OPTRAN2 and KFB for a Typical 
Section; M = 0.8, k = 0.01 
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M = 0.8, k = 0.001 
Pitch about leading edge 
KFB 

OPTRAN2; 72x60, |x| £ 10 
0 |y| ^ 10 

x jyj ^ 10 

□ |y| ^ io 



- 1.0 - 0.6 - 0.2 0.2 0.6 1.0 


x 


Figure 11. Comparison of Pressure Distributions From OPTRAN2 for Three Upper/Lower 
Boundary Locations and from KFB for a Typical Section; M = 0.8, k = 0.001, 
Imaginary Part 
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M = 0.8, k = 0.01 

Rectangular wing, aspect ratio = 3.0 
Pitch about leading edge 
Root chord 

— RH04 

OPTRAN3; 48 points in x, |x| ^ 8.0 
32 points in z, |z| ^ 8.0 
0 15 points in y, y ^ 7.0 
x 16 points in y, y 11.75 
□ 17 points in y, y ^ 16.75 

OPTRAN3: 

* 48x14x32; -8.0 ^ x ^ 3.0, |y| £ 4.5, |z| ^ 4.5 


Imaginary 


Figure 13a. Comparison of Pressure Distributions From OPTRAN3 for Four Outboard 
Boundary Locations and From RH04 for a Rectangular Wing; M = 0.8, 
k = 0.01, Root Chord 


45 





M = 0.8, k = 0.01 

Rectangular wing, aspect ratio = 3.0 
Pitch about leading edge 
r, = 0.885 

RH04 

OPTRAN3; 48 points in x, |x| ^ 8.0 
32 points in z, |z| <, 8.0 

0 15 points in y, y <, 7.0 

x 16 points in y, y <, 11.75 

□ 17 points in y, y <, 16.75 

OPTRAN3: 

# 48x14x32: -8.0 <, x <, 3.0, |y| <, 4.5, |z| ^ 4.5 


-5 -r~ 
- 1.0 



Figure 13b. Comparison of Pressure Distributions From OPTRAN3 for Four Outboard 
Boundary Locations and From RH04 for a Rectangular Wing; M = 0.8, 
k = 0.01, v = 0.885 




M = 0.8, k = 0.01 

Rectangular wing, aspect ratio = 3.0 
Pitch about leading edge 
Root Chord 

— RH04 

OPTRAN3; 48 points in x, |x| ^ 8.0 
32 points in z, jz| £ 8.0 
0 15 points in y, y ^ 7.0 

x 16 points in y, y ^ 11.75 

□ 17 points in y, y ^ 16.75 

OPTRAN3: 

# 48x14x32; -8.0 ^ x ^ 3.0, |y| £ 4.5, |z| £ 4.5 


imaginary 


-1 .24 — 
- 1.0 



Figure 14a. Comparison of Pressure Distributions From OPTRAN3 for Four Outboard 
Boundary Locations and From RH04 for a Rectangular Wing; M = 0.8, 
k = 0.01, Imaginary Part, Enlarged Scale, Root Chord 




Figure 14b. Comparison of Pressure Distributions From 0PTRAN3 for Four Outboard 
Boundary Locations and From RH04 for a Rectangular Wing; M = 0.8, 
k = 0.01, Imaginary Part, Enlarged Scale, rj = 0.885 
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M = 0.8, k = 0.001 
Pitch about leading edge 

OPTRAN2; 72x60, |x| £ 10 
+ |z| 10 

0 |z| ^ 50 

x |z| £ 100 


6 
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2 
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Velocity 
potential 
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Figure 15a. Comparison of Velocity Potential Along a Column of Points From OPTRAN3 for 
Three Upper/Lower Boundary Locations for a Typical Section; M - 0.8, 
k = 0.001 
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M = 0.8, k = 0.001 
Pitch about leading edge 


OPTRAN2; 72x60, |x| £ 10 

+ |z| £ 10 

0 |z| £ 50 

x |z| <, 100 



Figure 15b. Comparison of Velocity Potential Along a Column of Points From OPTRAN3 for 
Three Upper/Lower Boundary Locations for a Typical Section; M - 0.8, 
k = 0.001, Enlarged Scale 
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Figure 16. Steady Pressure Distribution Over the Root Chord of a Rectangular Wing With a 
NACA 64A010 Airfoil Section; Aspect Ratio = 3.0, M = 0.8 
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y 



Figure 1 7. Steady Pressure Distribution Over the Upper Surface of a Rectangular Wing With 
a NACA 64A010 Airfoil Section; Aspect Ratio = 3.0, M = 0.8 
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+ 

X 

0 


FAST 

0PTRAN3; with thickness 
0PTRAN3; without thickness 


M = 0.8 

Rectangular wing of Reference 1 9 
Aspect ratio = 3.0 
NACA 64A010 airfoil 



Figure 18c. The 0 2J Generalized Force Versus Reduced Frequency for a Rectangular Wing; 
M = 0.8, Aspect Ratio = 3.0 
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M = 0.8 

+ 

FAST 

Rectangular wing of Reference 1 9 

X 

OPTRAN3; with thickness 

Aspect ratio = 3.0 

0 

OPTRAN3; without thickness 

NACA 64A010 airfoil 






Imaginary 


Figure 18d. The Q 22 Generalized Force Versus Reduced Frequency for a Rectangular Wing; 
M = 0.8, Aspect Ratio = 3.0 















Figure 21. Comparisons Between Flutter Calculations From 0PTRAN3 and FAST; Flutter 
Frequency Ratio Versus Mass Ratio 
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APPENDIX 


A COORDINATE TRANSFORMATION SYSTEM FOR SWEPT AND 

TAPERED WINGS 


A.1 THE BASIC TRANSFORMATION EQUATIONS 

In setting up mesh patterns for swept and/or tapered wings, it is desirable to set up grid pat- 
terns with the points aligned in some sense with the leading and trailing edges of the planform. 
The formulation that follows provides such a mapping system starting with a rectangular region in 
the physical x-y plane. The planform leading and training edges are extended to the outer mesh 
boundary (fig. 22), and transformations between (x,y) of the physical plane and (£,??) of the transfor- 
mation plane are defined for the regions upstream of the leading edge, between the leading and 
trailing edges, and downstream of the trailing edge. The transformations provide for continuous 
second derivatives everywhere. Also, the sharp apex at the root of a swept wing is rounded to avoid 
singularities at the root plane. 

Let x = X LE (y) represent the leading edge of the wing, including its extension to the outer mesh 
boundary. The semichord S(y) is also extended to the outer mesh boundary. Both X LE (y) and S(y) 
have continuous second derivatives. For the wing region, - 1 < £ < 1 , we use the transformation 


£ = [*' X LE (y)]/S(y) - 1 

(A-l) 

v = y 

(A-2) 


Solving equation (A-l) for x yields 

* = X LE (y) + S(y)-(| + 1) 

differentiating with respect to x yields 

1 = SC VK, 


and 


= 1/SOO 


Differentiating equation (A-3) with respect to y leads to 

o = x' E + s'o-ws+D + i y n x 


(A-3) 


(A-4) 
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blank 


and 


£>■/£* — ~X LE — S'O0-(£ + l) (A-5) 

For the mapping upstream of the wing, $ < - 1 , we consider a form that has continuous second 
derivatives everywhere. Analogous to equation (A-3), we write 

x = X LE (y) + SOO-G+U + 3 /O') (A-6) 

We choose /O') so that at £ = £,, x = x,. Thus 
I *i = X LE + S(y)-($! + l) + /O') 

and we obtain 

x = X LE 0') + S(y) - (? + l) + x i • _ XleO') - SO') , ^$i + l) (A-7) 

We now find £ x and $ r Differentiating equation (A-7) with respect to x yields 
1 = ^S(y) + [*i - XleO') ~ S(y)-(s,+ l) 

and “ ‘/(SO) + [*i - X LE (y) - SOO‘tt + 1)]] (A-8) 

Differentiating equation (A-7) with respect to y yields 

o - + S'W-IS + l) - (l^) [x' E +S'w({, + l)] + {,/{, 

and l,/l, - -X^ - S'Wd+l) + (iii) [x^, + S-Ck)-({, + i)] (A-9) 

Downstream of the wing region where $ > 1 , we define a mapping of the form 

* = X LE + S O') • (1 + 1) + g(y) 

Ysmax 1/ 

and determine g(y) so that at $ = Z max we obtain x = x max . Then 

gOO = *max - X LE (y) - S(y) • (Imax + 1 ) 

and x = X LE 0') + s O') • (I + 1) + Y t ) ^max - XleOO - SO') • ($ ma * + l) (A-10) 

I 

I 
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Differentiating with respect to x yields 


i = soo • i, + * <s ■ _£ [*«.-x u O')-so>)- ({_+»)] i, 

vs max 

and $, = 1/ js(y) + ^ "Jjp [x max - X LE C v) - S(y) • (s ma , + 1)] J 
Differentiating with respect to y leads to 

0 = X' E + S'tvMS + l) - (^rr) 3 [ X LE + s, Cv)-(«max+l)] + €,/€* 

and %y/% x - -X^ E — S'0')*(? + 1) + ^ ^XleCv) - S'(y)-^ max + l) 

In the finite difference operators, we need 

f = i n x , G = Z y / S x 

and 

g 2 = g 2 /f 


(A-ll) 


(A- 12) 


which are defined by equations (A- 4) to (A-12) above. The metrics are calculated from these 
expressions in the program. 
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A.2 EXTENSIONS OF THE LEADING AND TRAILING EDGES OF THE WING WITH 
SECOND-DERIVATIVE CONTINUITY 


To eliminate the singularity of the leading edge at the line of symmetry, y = 0, we turn the 
leading edge so that it has zero slope at the plane of symmetry. Thus we express the leading edge in 
the form 


X L e = X le (2) + af + by' for y(2)<y<y(N) 

(or 0<y<y N ) 


(A-13) 


where y(2) = 0 is the plane of symmetry and y - y(N) as y N is the value of y at which the slope of 
the new curved edge matches the slope of the straight leading edge (see fig. 22). Thus we require 
X leO'n) = 0, since we consider wings with straight edges. Thus 

Xle = 2 a + 6 by^ = 0 


and 

b = -a/(3y N ) 

Then X LE (y) = X LE (2) + a[y 2 -y } /( 3^)] 

At y = y N , we have X^ £ = y N tan 0 LE = c,. Differentiating equation A-15 then yields 


(A- 14) 
(A-15) 


X le(/n) = «[2>' n -3>'n/(3> , n)] = 

Hence 

a = c,/y N 

The equation for the leading edge for the section near the plane of symmetry then becomes 

X LE = X le (2) + c , [y 2 / y^—y 3 / (3.Vn)] for 0<y<y N (A-16) 

At y = y N , we have from equation (B-49) of Reference 8 


XleCvn) = - 1 + 


Thus equation (A-16) at y = y N becomes 

X L e = X LE (2) + 2c,y N /3 
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m boundary 



i Leading and Trailing Edges of a Swept and 
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Hence a = c,/y N and 


XleCVn) — X le (2) + 2c t y N /3 — — 1+c^n 
Finally, we obtain for equations (A-16) 

X LE (y) = -1 + c,[y/y N -y 3 /(3yN) + Jn/ 3] forO<y<y N (A-17) 

Similarly, for the trailing edge we write 

X TE (y) = 1 + (c, + 2c 2 ) • [/An -//(3y N ) + y N /3] (A-18) 

from which we obtain 

SCO = [x TE (y) - X LE (y)] / 2 = 1 +c 2 [y 2 /y N -y/(3y N ) + yjs] (A-19) 

For y, >y>y N , i.e., the straight portions of the leading and trailing edges, we have 
from equations (B-49) and (B-50) of Reference 8 

X LE O0 = -l+c,y (A-20) 

X te O) + 1 + (q + 2 c 2 )y (A-21) 


Extending the leading and trailing edges beyond the tip, we find that the two edges intersect 
at y=y, 

or 

-1+c,^ = 1 + (c, + 2c 2 )y t 
or 

y, = - 1 /c 2 (A-22) 


Differentiating equations (A-17) and (A-14) with respect to y yields 

X LE O0 = c, (2 t-t 2 ) , t = y/y N (A-23) 


S 'O') = c 2 (2t-f) 


(A-24) 
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For y N < y < y t t we have 


XleCv) = -i + qy 

(A-25) 

XteCv) ~ 1 + (o, + 2c 2 jy 

(A-26) 

S(y) = 1 +c,j> 

(A-27) 

Xle O') ~ c i 

(A-28) 

S ' O') = c 2 

(A-29) 


We extend the wing leading and trailing edges beyond the tip, y„ with a cubic in y. Choose 
the point y = y m to lie halfway between wing tip and the intersection of the linearly extended 
leading and trailing edges, or 

y m = (?,+y t )/2 - {y-\/c z )/2 (A-30) 

Thus, for y, < y < y m , we write for the equation of the leading edge 

XleOO = -\+C\y+b{y-y,f (A-31) 

The mean slope is chosen for both leading and trailing edges at y = y m . Thus 

[(c|+2c 2 ) + 0,1/2 = 0(4-02 

Hence at y = y m 

X LE (ym) = c x + 3b (y m -y ! f = c, + c 2 

or b - cj [3(y m -y,) 2 ] 

Similarly, for the extension of the trailing edge we write 

XteOO = 1 + (c, + 2 c 2 ) y + by (y-y,) 

XteCf) = (o, + 2c 2 ) + 3 by(y m -y,f = c,+c 2 
Thus, at y = y m , by becomes 

b\ = -c z /[i{y m -y,f] = -b 


(A-32) 


(A-33) 


(A-34) 


(A-35) 
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After substituting b and b u the leading- and trailing-edge extensions then become for 


y,< y < y m 

X le O) = - 1 + cj + c 2 (y -y,f /[3(y m -y,f ] (A- 36 ) 

X te O) = 1 + (c, + 2 c 2 )y - c 2 (y-y,f/\ 3 (y m -y,) 2 ] (A-37) 

For the semichord SO) we obtain 

SO) = l+c 2 {y-(y-y t f/\3(y m -y,f ]) (A-38) 

Differentiating equations (A-36) and (A-38) with respect to y yields 

X le O) = c, + c 2 (y -y,) 2 / (y m -y,) z (A-39) 

SO) = c 2 [i-(y~y t y/(y m -y,) 1 } = c,+c 2 -X te O) (A-40) 


To simplify the outer mesh boundary conditions, we extend the leading and trailing edges 
for y = y,„ to y e , the outside boundary, by a third-degree polynomial, so that slopes 
X' = X' = S' = 0 at y e . Thus we write 

LE TE 

XleO) = x LE m + a \{y-y e f - (y m -y,f ] 

, f / >3 r ,3 ] (A-41) 

+ b [ {y-y e f - \y m -y t ) | 


where X LE m is the x leading edge coordinate at y = y,„, the point halfway between y, and the 
point of intersections of the leading and trailing edges. From equation (A-36) we have 


X le O) = -1 + c,y m + c 2 (y m -y,)/3 + a [(y-y,) 2 -(y m -y e ) 2 ] 

+ b [(y-y e ) 3 - {y m -ye) 3 } 


(A-42) 


At y, = y m , we must have X^ E = c, + c 2 , or the average of the leading- and trailing-edge slopes at 
the wing tips. Thus 


XusOm) = a[2 (y m -y e ) + b[3(y m -y e f] = c, + c 2 (A-43) 

To match the second derivative at y = y m , we have from equations (A-39) 

XteOJ = 2 c 2 /(y m -y) 

Thus from equation (A-42) we obtain 

2a + 6b(y m -y e ) = 2 c 2 /(y m -y,) 

or 

a + 3b(y m -y e ) = c 2 /(y m -y,) (A-44) 
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Eliminating b between equations (A-43) and (A-44) yields for a 

a = (c, + c 2 )/ (y m -y e ) - c 2 /(y m -y,) 

From equation (A-44) we have 

3 b(y m -y e ) = c 2 /(y m -y,)-a 

or b = \c 2 /(y m -y,)-a\/\3(y m -y e )\ 

Similiarly, for the trailing edge we write 

X TE (y) = 1 + (c, + 2 c 2 )y m - c 2 (y m -y,)/3 

\{y-y e f - (. y,r,-y e f | 

+ b x \(y-y e f-(y m -yef | 

At y = y,„, X' XE = c, + c 2 and equation (A-47) yields 

2a, + 3b x (y m -y e ) = (c, +c 2 )/(y m -y e ) 

At y = y„„ we have from equation (A-37) 

XteOO = ~ 2c 2 / (y m -y,) 

Thus 

X TE (y) = 2a { + 6b x (y m -y e ) = -2 c 2 /(y m -y,) 

Oi + 36, (y m -y e ) = -c 2 /(y m -yi) 
Solving equations (A-44) and (A-45) for a, and b t yields 

= (C| +c 2 )/(y m -y,) + c 2 /(y m -y,) 

and 

b i = -\c 2 /(y m -y,) + a l \/\3(y m -y e )] 

Let 

X T em = 1 + (Ci + 2 c 2 )y m - c 2 (y m -y,)/ 3 
-ai(y,„-y*) 2 - 6 1 (y m -y ? ) 3 


Then 


X TE (y) - X TEM + a\(y-y e f + b x {y-y e f 


(A-45) 

(A-46) 

(A-47) 

(A-48) 

(A-49) 

(A-50) 

(A-51) 

(A-52) 

(A-53) 
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Let 

X L em = -1+c \y m + c 2 (y m -y,)/3 

®(.Vm ye)‘ ~ ^[y hi y ? ) 


then 


XleOO = X LEM + a (y-y e f + b{y-y e f 


(A-54) 


(A-55) 


Let 

I 


Slem - (Xtem - Xlem)/2 


(A-56) 


then 

S(^) = S LEM + (a^-a^y-y,) 2 / 2 + (ft, -b)(y-y e ) l /2 

X LE Cv) = 2a(y-y e ) + 3b(y-y e ) 2 (A . 57) 

S'OO = (a,-a)(y-y e ) + 3(b,-b)(y-y e ) 2 /2 

For convenience, we summarize the results and include the first derivatives of the leading 
edge and semichord that we will need in the calculation of the coefficients of the difference 
equations. 

For 0 <y< y N 

X LE (y) = -1 + c,y^(t 2 ~t } /3 + 1/3) , t=y/y N (A-58) 


S(y) = 1 + c 1 y ii (r-t i /3 + \/3) 

(A-59) 

X LE Cv) = c,(2t- 1 2 ) 

(A-60) 

S' O') = c : (2t-t 2 ) 

(A-61) 
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For y^ < y < y, 


XleCv) — ~ i + c i y 

(A-62a) 

X TE Cv) = SCO + (cj + c 2 )j 

(A-62b) 

S Cv) = 1 + c^y 

(A-63) 

XleCv) = c i 

(A-64) 

S'Cv) = c 2 

(A-65) 

For y, < y < y m 


X LE Cv) = -1 + c l y + c 2 (y-y,) 2 /\3(y m -y l f | 

(A-66) 

S (y) = 1 + c 2 (y-y,) i /\My m -y l ) 1 | 

(A-67) 

X LE Cv) = C\ + c 2 {y-y t ) 2 /(ym-y<f 

(A-68) 

S'CV) = c l+ c 2 — Xle 

(A-69) 

For y m < y < y e , where y e = .y(JMAX) 

X L em = - 1 + c,y„ + c 2 (y m -y,)/ 3 
-a(y m -y e f-b(y m -y e ) 3 

(A-70) 

XleCv) = X lem + a(y - y t f + b(y - y,)* 

(A-71) 

a = (c,+c 2 )/(y m -y e ) - c 2 /(y m -y,) 

(A-72) 

b = [c 2 /(y m -y,)-a]/[3(y m -y e )] 

(A-73) 

X TE m = 1 + (C] + 2c 2 jy m — c 2 (y m — _y,)/3 

(A-74) 

-My* -ye? - bi(y m -y e f 


XteCv) = X TEM + a\[y—y e f + ^i(.v — .v*) 3 

(A-75) 
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a, = {ci + c 2 )/(y m -y e ) + C 2 /(y m -y,) (A-76) 

by = - [c 2 /(>' m ->',) + a 1 ]/[3(>' m ->' f )] (A-77) 

Slem = (^tem - -^lem)/2 (A-78) 

SOO = Slem + (a, -a){y-y e f/2 + (by- b)(y-y e f /2 (A-79) 

X LE Cy) = 2a(y-y e ) + 3b(y-y e f (A-80) 

S 'O') = (a l -a)(y-y e ) + 3(b l -b)(y-y e j 1 /2 (A-81) 
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A.3 TRANSFORMATION OF THE DIFFERENTIAL EQUATION 
FOR A SWEPT AND TAPERED WING IN TRANSONIC FLOW 

A.3.1 Basic Formulation 

The governing equation for the unsteady potential for harmonic motion of a wing in transonic 
flow is given by 


(u<P\x)x + <P\yy + <P\zz~ 2iu, filx/t + 0>V|A = 0 


(A-82) 


where u = K -(7+ \)<p ox . The quantities <p 0 and <P\ are the steady and unsteady perturbation 
potentials, respectively. The variables x, y, z are the Cartesian coordinates. We define a coordi- 
nate transformation that aligns the leading and trailing edges with coordinate lines. For 
points on the wing, we let 


£ = [x- XleOO] / S(y) - 1 


n=y 


(A-83) 


r=z 


where X LE (y) is the x coordinate of the leading edge as a function of the spanwise variable y 
and S(y) is the semichord. We consider a more general transformation 




v =y 


(A-84) 


t=z 


then 


<Plx = ‘Pl^x ‘Ply = ‘Plr, + 'Plfiy 

The differential equation becomes 

£x( u <Pl(£x)( + (^ 1 , + + £y(‘Pir, + < P\£y)l 


(A-85) 


(A-86) 


+ v J irf-2/wvJ U ^A + ojViA = 0 
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To express the equation in conservation form, we consider 



(A-87) 


Since x = x(£»jO differentiation with respect to x and y yields 

1 = =x£ x or £r=l/x* (A-88a) 


0 = X{£> + X, 


Hence 


or 


I y = “X n /X| 



(A-88b) 


(A-89) 


When the two equations (A-88) and (A-89) are combined in equation (A-86), the first terms on 
the right-hand side of' equation (A-87) cancel and equation (A-86) becomes 


+ <®irr/£* “ 2 iuxpi ( /e + «Vi/(e£*) = 0 (A-90) 


Let F = 1 /£, 0 = Z y /Z x d 2 = t)/$ x w = w|. 
Then the differential equation takes the form 


[(u + G 2 )<p u ] { + (F<pi,), + (F<p, r ) r 

+ (G<p,,)j + (Gy> u ), - 2io)<fi u / e + u 2 Fip l / 6 = 0 


(A-91) 
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A.3.2 The Derivation of the Difference Equations 


We now express the differential equation in difference form, following the notation of Reference 
7 for the coefficients. For the first term, we have for elliptic points 

2 = C Mi+\/ljk {Vi+\jk~ Vijk) ~ di^i- 1/2/* {’Pijk ~ ’Pi- Ijk) (A-92) 

and for hyperbolic points 

2 = C <- l“i- W 2 jk {f ijk ~ Vi- ijk) ~di-l U i-}/ 2 jk {Vi-\jk~ Vi-ljk) (A-93a) 


where 

c,- « l/[(*j +1 - *,_ ,)(*,+ 1 - *,)] 

and (A-93b) 

d t = l/[(x i+l - *,_,)(*, - *,_,)] 

Following the procedure in the ADI method, we combine the two procedures and obtain 


D2l 


{Vi+\jk V ijk ) \/2jk {Vijk Vi- Ijk) 


(A-94a) 


+ D1X*. [c,-_ x Uj_ W2j k (v ijk “ Vi- \jk) ~~di- (Vi-\jk~ Vi-ljk) 


The variables D2 and DlX are zero or one for elliptic or hyperbolic points to select the appropriate 
operator. For a shock passing between points ( i,j,k ) and (i + l,j,k); 


(see ref. 7). 

For the second term, we write 


DlX* = (x, - x i _ 2 )/[(x i - x i _ 2 ) + (x i+l - *,_,)] 
D2 * = (x /+l - Jf,_ ,)/[(*■, - x,_ 2 ) + (x i+l - *,_,)] 


(A-94b) 




— c fi*2i+WZjk (jPi+\jk~ < Pijkj~ djQn-\/2jk (Vijk ~ Vi- Ijkj 


(A-95) 
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For the remaining terms, we have 


2 £®V,j — Oyj ^'ij-\/2k{fij-\k~'Pijk) t>yj^ij+W2k {fijk ‘Pij+lk) 

\ (*Vr)r = F .y [<*Jv>,y*- »-*>«/*) " b * (*«» “*<»+>) 

^ (G^i,) = c iA+iy J^ c i///+iy+u + (^|>> -C i>y) 'Pi+ijk-diyjVi+ij- 
+ (d\, _ C l:j Gy ^Cj^y + 1* + (d\yj ~ C\yj) <Pjjk ~ d\yjPi- \j- 1*] 

— d u Gj_ jy ^i^j-iy+1* + {d\yj— Cly^Vi-ijlc — d\y<Pi-ijk 

— [<*,] = C\y- Gy + 1 j^ C U<^,+ 1_/ + 1* + ¥ J i/+l* — d\yjir i- 

+ {d X yj C\yj) G,yj^C| ( ^ (+ ijic + (d u - c nj <p uk - d u <Pj_ i jlc 

~~ d\yj Gy- 1 \j- 1 *• + {d\i~ C\ij <Pij-\k~ d\j<Pj-\j 


/-iy- 1 * 


For the first derivative with respect to £, we have for elliptic points 

- iuxp | ^ / e = - (/«/ e) [c u (<Pi+ ijk ~ <Puk) + d\i (‘fiijk ~ ‘Pi- 1 jk) 
and for hyperbolic points 

-i(*Pn A = “ (*«/ e) [c a (?</* -V»/-iy*)-^i/-i (<Pi-ijk~ Vi-iik) 
At the shock point, we use the operator defined in Reference 7. Thus 



<Pi + \ik + ‘Pilk - 'Pi- I ik - 'Pi-lik 

yu/e) 

£i +1 + £i - £i- 1 - £(-2 


We combine these three operators in the form 

— iuip^/e = (/ 0 >/e) \jk + ^hkl'Pijk + C hki'Pl- ijk + C hk‘\'Pi- 2 jk 


(A-96) 

(A-97) 

(A-98) 

(A-99) 


(A-100) 
(A-101) 
(A- 102) 
(A- 103) 
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where the coefficients are listed in the table below 


C /!*l 

c hk 2 

C A*3 

C / 1*4 

Elliptic -C| 

l l 

*>i 

0 

Shock point - 1/A£ 

— 1/A£ 

+ 1/A£ 

+ 1/A£ 

Hyperbolic 0 

“3 2, 

c 2/ + ^l/-l 

- ^li-I 

where 





= S/+ 1 + S< ~ l/-i ~ Si-2 


Finally, combining equations (A-94) through (A-103), we obtain the following difference equation 
D1X* ^C,-_ |W,_ X /Tj k (<Pijk ~ 'Pi- \jk) ~ di- l^i-3/2 {pi-\jk~ Pi-Xik 

■b D2 k 1 /2yit (fPi+\jk~ 'Pijk)~ diUj-\/ljk (pijk~ 'Pt-Xjkj 

+ C @2i+ 1/?/ (V,+ 1 ;* - ^i/Vt) - dftli- Vlj ('Pijk — 'Pi - ly*) 


“b y_ 1/2 (pij- \k 'Pijk) byjF ij+ 1/2 (pijk < Pij+\k) 


+ F„ 


■iy ^V*) (^<7* 'Pijk+l) 

"b |^C| ( G I+ |y j^ly^i+ly'+l* "b — C !>’_,■) 'Pi+\jk~ diy/Pi+lj-lk 

"b (^J( _ C|/) G,y + I* + {d\y.— C\y)jiPijk ~ d\y.<Plj_ 




l ] /2 


C\yjPi- 1 7 + 1 * "b (d\yj c lyj) 'Pi- 1 jk d \yjPi - \j - 

■b^lyy^(/+t [ c l(^/+ly'+l* "b (dy — ^-u)Pij+ Ur — ^\iPi- ly + Mr) 

+ (^l>y -c lyy) G/y [ c l/V’/+lyit + (^1/ — c l/)^y* - ^l/V’i- ly* 

— d\ y Gjj_ 1 J^c 1( v?; + jy_ n- + (dy— Ciij'Pij-ik — d\i<Pi-ij~\k ^ /2 
+ (/w/e) [c A/tl <p l + lyVt + C hk2 <fiijk + Chki'Pi- \jk + c hk*'Pi- 

+ w 2 F ij <p i j k / (2e) = 0 


(A-104) 
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In figure 1, we illustrate the labeling of the coefficients of equation (A-104). With this labeling, the 
difference equation takes the form 

SUB(*)*?0fc-i + DIAG(Ar)*^ + SUPER(Ar)*v>,y* + 1 + CFOR2(A:)^_ ?/ * 

+ CFOR(k)*<p i _ ljk + CAFT(k)*<p i+ ijk + CFORIN (*W,- v -u (A-105) 

+ CIN(k)*^. u + CAFTIN(Xr)v, + 1 y-u + CFORO(Ar)V,- V+I * 

+ COUT(it)*^ +1 * + CAFTO (k)*v i+vk = RHS(Ar) 

Combining all like terms in equation (A-103), we obtain for the coefficients 

'Pi- ij-ik' d u d\yj (G,_jy + G,y_j)/2 = CFORIN(Ar) 

Pij- i* : a y^ij-\n~ d\ y j (</„•- c u ) (Gy Gy_i)/2 = CIN(Ar) 

Pi+ ij- 1*: -Cudiyj (G, + iy + G y _ 1 )/2 = CAFTIN(k) 

<Pi- \jk* D2^,« 

/- 1/^/Vt ” D 1 X, (Cj_ x u t _ l/ y k + d,_ iUj-w2jk) + ,a)C A«/ e 

+ c/,G, ( . -]/2j — d\ id\y. — c i>y)(G,'- iy + Gyj/2 * CFOR(Ar) 
v^y*' — D2* (c i Uj +l/2 j k + djiij. i/yk) + + ( ,w / e) c AA . 2 

ij- i/2 ~~ b y F tJ+ 1/2 — F jj (Ozii "b 
— c f*7i+l/7j — d£x‘2i-\/y + (^li — c l;)(^l>y — c ly y -)Gy 

+ w 2 F,y/ (2e) = DIAG(Ar) 


*Pi+ 1 jk* 

^2k C i^i+ 1/2/* + c @2i+l/2j + J ( w / e ) c hk\ 

+ C U (di y .-Ci yj ) (G i+ ij + Gy)/2 = CAFT(Ar) 

(A- 106) 

*Pi- \j+ \k' 

— d\iC\yj (Gy + 1 + G,_ iy)/2 = CFORO(Ar) 


<Pij+ \k : 

byj^ij+ 1/2 + C l>y (<*W — C li) (Gy + Gy + ( ) /2 = COUT(At) 


<Pi-2jk' 

BlX k d j _ l u i -^ /2jk + i(j)C hki / e = CFOR2(Ar) 


*Pi + \j + \k • 

c \yf\i (G/ +v + G tf+1 )/2 = CAFTO(Ar) 


A.3.3 Boundary Conditions on the Mesh Boundary 


On the left mesh boundary, we apply the outgoing wave boundary conditions used in References 3 to 

8. Thus, on x- 

(x, + x 2 ) / 2, we have 



- /wlVW(l-M) = 0 (A-107) 
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7f 


SUPER(K) 


i, j, k + 1 



Figure 23. Definition of the Coefficients for the Finite Difference Operator 


Since % x = 1/F, we obtain for the difference form of equation (A-107) 


tnkZJjm. _ k»MF 3/2/ <Pm + V\ik _ o 

1-M 2 


or 


1- 


/cjM-A^rF 3/2 ,- 


2(1 -M) 


<P2jk 


' + 2(1 -M) J*"'* 


= 0 


where A£, = £ 2 ~ £i We write 


where 


and 


^iy*r - C^jVljk 


^K1 j - 


1 ~ C KI/ 

1 + C Kl j 


_ /oiM-ASi „ 

Cku " 2(1 -M) ),v 


Similarly, on the downstream boundary we have 


<PiiZx + 


1+M 


• <p = 0 


In difference form, this becomes for i = i mttX — 1 


Vj+nk 'Em j_ £^M£l± i Hi . < £i+Aj!L^JBiik — q 

£,+ 1 ~%i 1+M 2 


or 



/coM- A£ m *F /+1/2 y 


/cjM * A£ m * F, + 1 /2 / 

*Pi+ 1 jk 

1 + 2(1+M) J 

Pijk 

2(1 +M) 


Introducing similar notation leads to 


^K2/^/ miut \jk 


(A- 108) 


(A- 109) 


(A-110) 


(A-lll) 


(A-112) 
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where 

^K2J = (l “ C K2>)/(^ +C K 2j) (A-113) 

c K2J = /coM-A£ m •F w _, /w /[2(1 + M)] 

On the line of symmetry, y = n = 0, we apply the symmetry condition 

<Piy = <fiu • 0 

In difference form, this becomes 

v^/i* = ^(3* (A-114) 

where j = 2 is the plane of symmetry and leading and trailing edges are turned normal to the plane 

y - X2) = 0. 

At the outboard spanwise mesh boundary, y = + >y max -i)/ 2> we apply the outgoing wave-type 

boundary conditions in the form 


/«mVk 


+ T^M r *’ 1 = 0 


or 


„ /cjMVk 

¥> i , + €^ 1 | + J _ M 2 ^1 = 0 

We select a mapping so that £ „ = 0. Hence, for j = y max - 1 , we have for the difference equation 


This may be written as 


where 


< £ll±AK *Piik 

Vj+\ ~~ Vj 


/ojMVK \ <g M±lk ±<p M 

1 -M 2 2 


= 0 


'Pij+ l* - C K3 Vijk 


C K3 ~ (1 -C K3)/(1 +C K3) 


c K3 = ZwMVK-A^/ 2(1 -M ; 

= Vj maJC ~ Vj max - 1 


(A-115) 


(A- 116) 
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On the lower boundary where z - (z t + z 2 )/l, we have 




iu MVK 

= 0 


which leads to 

‘Pijl — ^K4 <Pijl 

where 

^K4 = (1 ~ C K *}/ (1 + C K4) 

c K4 = /«M • Az, -Vk / |^2(l -M 2 ') 

. Az, = z 2 -Si 


(A-117) 


(A-118) 


Similarly, on the upper boundary where z - (z k _i + z k )/ 2, we have 


+ 


toMVK 

1-M 2 


v? = 0 


or 

^V*max = ^K*/^y* max - 1 (A- 11 9) 

where 


^K 5 (1 "“ C Ks)/( C K 5 ) 


c K5 = • VK/[2(1-M 2 ) 


^ Zm Zk max Zk max ~ 1 


(A- 120) 


We now apply the boundary conditions on the mesh boundaries to the difference equations. On the 
left boundary * = (*,+ x 2 )/ 2, the term 


CFORIN(Ar) • <p {j _ lk 


is replaced by 


CFORIN(*)- *>*,._,* 
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then these two coefficients are modified by 


CIN(Ar) = CIN(Ar) + c^-CFORINfA:) 
CFORIN(it) = 0 

Similarly, CFOR (k) ■ ip l]k is replaced by CFOR (k) • c Ki ^ 2j k 
and the coefficients are modified by 

DIAG(Ar) = DIAG(Ar) + c Klj * CFOR(Ar) 

CFOR(£) = 0 


On the plane of symmetry we have, for y = 2, 


< Pi\k ~ ilk 

and the coefficients are modified by 

COUT(Ar) = COUT(Ar) + CIN(Ar) 

CIN(k) = 0 

On the boundary outboard of the wing, we have for y = y max - 1 

< Pij+ 1* = C K 3 <Pjjk 

and the coefficients are modified by 

DIAG(Ar) = DIAG(£) + c K3 * COUT(Ar) 
COUT(Ar) = 0 

Similarly, on the right, downstream boundary, we substitute 

^W* = C K?/^W-l (/* 


and modify the coefficients by 


DIAG(Ar) = DIAG(Ar) + c K2j * CAFT(Ar) 
CAFT(fc) = 0 

CIN(Ar) = CIN(A:) + c K2j * CAFTIN(Ar), y * 2 


(A-121) 


(A-122) 


(A-123) 


(A-124) 


(A-125) 


(A-126) 


(A-127) 
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CAFTIN(Ar) = 0 


(A-128) 


COUT(Ar) = COUT(Ar) + c K2j * CAFTO (*), j*j max - 1 

CAFTO(Ar) = 0 

On the lower boundary for z - (Z\ + z 2 )/l and k =2, we have similar conditions, namely, 

V’i/I = ?K4 <Pij2 

and the coefficients are modified by 

DIAG(2) = DIAG(2) + c K4 *SUB(2) (A-129) 

SUB(2) = 0 

Similarly, for the upper boundary we obtain 

^*max = C K5 ^max " 1 (A-130) 

and 

DIAG(Ar max - 1) = DIAG(A: max - 1) + c KJ * SUPER(/ max - 0 (A-131) 

SUPER(Ar max -l) = 0 
A.3.4 Wake Boundary Conditions 

The condition that the pressure be continuous across the wake is given by 

A (p x + iuAifi = 0 (A- 132) 

where A<p is the jump in unsteady potential across the wake. This may be integrated to yield 

A?, = Atp i{ + , exp( - iu(x,-x, x + ,)) (A- 133) 

We express the argument 

ARG = -a> (*,-*,• I + 1 ) (A- 134) 
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in terms of the variable £. Thus, for the downstream region for £ > 1, we have from equation (A-7) 
in the section on the coordinate transformation 


= X LE 0) + St/H€ + 1) + (—-37) 3 ^max - X LE (y) - S(J ) • (| max + 1) 

\smax A / L 


(A-135a) 


where we have chosen X max = £ max . For the trailing-edge point /, + 1 from equation (A-10), we obtain 


X i\ + 1 ~ X L e(T) + S(y) • ({,, + 1 + l) + ' ' + I ™' ) lmax - X LE 0)-S0')-(| max + l) 

\smax 1 / 


(A-135b) 


and for the general point x 


x i - X LE (y) + S(y) • (£,■ + l) + - \ $ ma x ~ X LE (y) - SO> (£ max + 1) 

\smax A / L 


(A-136) 


Thus the argument in equation (A-134) finally takes the form 

ARG - « [so>(£,, - X LE 0) - SO> (£„„ + 1)] [ (f^y - (j^-rr)’ ]] <A ' 137 > 

To satisfy the Kutta conditions at £ = we set the pressure jump at the trailing edge equal to zero, 
i.e., from equation (A-132), 

A<p lx + iu’Aifii = 0 


which in the transformed grid variables becomes 

A<fii( + iuF , A<fi i = 0 


(A- 138) 


At £ = or x = Xj , this becomes in difference form 


Ci/,-(A ^ 1 + 1> -A^ v ) + d U{ • (A<p ix j - A + iuiF^j- Aifi^j = 0 

We solve this expression for the value of A*?,, + , at the first point in the wake and obtain 
A *V'> = A *<V ~ ( rf i.,Ai«,)( A ^i> “ ( F <v/ C i/,)' A ^V 
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or 

^P/j + ly = (1 — ^\i\/^\i\ ~ i\J / ^\i\) (^lij/^I/'i) ^Piy — lj 


This may be written more conveniently as 

A<P/j + jy = c Kcl ’ + CKc2 * - ly 


where 


^Kct — 1 ^\i\/ c \i\ /wF /jy'/ci/j 


c Kc2 ~ ^Wj/ C l/[ 


(A-139) 


(A-140) 


Av5,y — *Pijk m + 1 i Pijk m C Sl (vijk m + 2 <fiijk m + [) 

~ C S2 ( — ^iy* m ~ ^ijk m - 1 ) ~ (^SlFy + ^S2^'(/) 

In general, the jump in potential across the airfoil is seen from Reference 1 to be given by 

&<Pij = <fiijk m + 1 - <fijk m ~ (-Slip ijk m + 2 ~ *Pijk m + l) 


i 

~ C S2 (V|y* m _ { Pijk m - l) ~ (^SlF,y + 


(A-141) 


where 

c sx = 1/[4S 1 (S 1 + 1)] c S2 = 1/[4S 2 (S 2 + 1)] 

^si = A(2S X + 1)/^4(S X + l)j d§2 = /?(2S 2 + 1)/[4(S 2 + 1)] 
h = ^ m + I - Z* m 

S 1 = ( Z *m + 2 _Z *«+») S 2 = { Zk m~ Z *m-') 
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Combining like terms leads to 


A<Pij = (l +Csi)'¥’y* m +l 


(l + c S2 ) • <p ijkm 

+ C S l<Pijk m + 2 + c S2'Pijk m ~\ 


(dntf + daFi) 


Then, to satisfy the Kutta condition, we have 


A^/,+ ly - C Kc l 


(l+Csi) Viykm+l ~ (l +C S2 ) Vi\jk m 

-C S1 Viyjkm + 2 + C S2 ( Pi l jk m - 1] 


+ c Kc2 


(l + C S1 ) </>/, _ ijk m -l (l + C S2) ^/'i - 1 jk m 
-C S l^/,-ly* m + 2 + C S2 <p ix _\jk m + ^ 

~ c Kcl ' (^SlF y + d S2 F y) “ C Kc2 ‘ (^Sl^ ^ - \j + ^2^ - Ij) 


For convenience, we define new variables 

All = C K< . 2 C S 2 

Ai 2 = -CKc2 (1 + C S2) 

\ 

A13 = C Kc2 (1 + C S l) 

A 14 = -Crc2 C S1 

A 2 i = c Kcl c S2 

A 22 = -^Kcl (! + C S2) 

A 23 = c Kcl (1 + c S i) 

A 24 = -CrciCsi 

Thus we have for 

A y>| 1+v = A u <^i,-ijfc m -i + A 12 <Pi x ~\jk m + A 13 

+ A 14 SO/, -ljk m + 2 + A 2 i >Pi i jk m -l + A 22 Vi\jk m 
+ A 23 <Ptyk m +l + A ■24<Pi l jk m *2 ~ c Kcl (^Sl F/p + d S 2^i u ) 
_c Kc2(^Sl^-l7 + ds2^'i l - ly) 


(A- 142) 


(A-143) 


(A-144) . 


(A- 145) 
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For the value of k below the wing surface k = k m , we add for the ith column for />/, + 

^ z k : m F yA^ij + U ®^P (ju(Xj — Xj^ + ] 

and define 

Wrii = - b z A U F ij 

m 

Wri2 = ~b Zk A 12 F ij 

ffl 

W K1 3 = _ b Zk A 13 F jj 
Wru = -b Zk A 14 F jj 

rn 

Wr2i = ~b z A 21 F ^ 

m 

Wr22 = ~ b z A^Fy 

K m 

Wr23 = -b Zk A 23 F ij 

n, m 

Wr24 ~ -b Zk A 24 F ij 

m 

RHSL = b Zkm j^c Kcl ^rf sl Fy + d a2 F yj + c Kc2 (d si F , t _ + d S2 F^_ Iy ) 
Similarly for k = k m + 1 we add for />/, + 1 and j less than the tip value 

♦ U eX P ,, + ,)) 

and define the variables 

W K31 = °*t m+ l A H F V 
W K32 = a ** m+ , A t2 F tf 
Wr33 = ^~k m + 1 A 13 F i>' 

W K34 = a z ,A 14 F /j 

W K41 = ° z k m + i A 21 F <>‘ 

Wr 42 = ^ m+1 A 22 F,y 

Wr4 3 = a z km+ 1 A 23 F (, 

W K44 = a *k m + 1 A 24 F V 

RHSU = -a z . + .-RHSL/ b z . 

K m K m 


(A-146) 


(A- 147) 


(A- 148) 


(A-149) 
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(A-150) 


Finally, we set E u = —exp ^ — — .*/_ 1 )) and for k — k m , i>i\ + 1 

RHS(KM) = RHSL • E u 
WAK11 = W KU • E u 
WAK12 = W K12 • E u 
WAK13 = W K13 • E u 
WAK14 = W KM • E u 
WAK21 = W K21 • E u 
WAK22 = W K22 • E w 
WAK23 = W K23 • E a 
WAK24 = W K24 • E u 
Similarly, for k — k m + 1 and »>i'i + 1» 

RHS(KMP) = RHSU • E„ 

WAK31 = W K31 • E m 
WAK 32 = W K32 • E u 
WAK33 = W K33 • E u 
WAK34 = W K34 • E u 
WAK41 = W K41 • E u 
WAK42 = W K42 • E u 

(A-151) 

WAK43 = W K43 • E u 
WAK44 = W K44 • E m 
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For k = k , the coefficients of $ associated with the wake are 

tn 

-6 -a*. = WAKll-^.-^ m -, + WAKl2^ tI _ 1An 

m 

WAK13 • <Pj x _ \jk m + 1 + WAK14 +2 

WAK21^ lV * m _, + WAK22 
WAK23 ■ <fii x jk m + 1 + WAK24-^, v * m + 2 

Similarly for k - k m + 1, i>h + 1 

a z , -AVij = WAK31 • _ y* m - 1 + WAK32 • <*, . 

+ WAK33 \jk m + i + WAK34- V 5, 1 _ uA r m + 2 

+ WAK41 • <Pi l jk m - 1 + WAK42 • <f>i x jk m 
+ WAK43 • ?, v * m + 1 + WAK44 • ^ + 2 


(A-152) 


(A-153) 
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A. 4 THE COEFFICIENT MATRIX AND SOLUTION PROCEDURE 


The resulting set of finite difference equations may be written in the matrix form as Ax = B. 
The coefficient matrix, A, for a rectangular planform is a sparse asymmetric matrix with 
coefficients arranged along the principal diagonal, four bands of unit width above the diagonal, 
and five bands of unit width below the diagonal. The number of nonzero elements in any row is 
nine for a rectangular planform. The actual bandwidth of the system is large, being 
NS*(JMAX-2)*KXB + 2*KXB, where NS is 2 for subsonic and 3 for mixed flow, and KXB is 
KMAX-2 for a full asymmetric flow and KM AX/2-1 for symmetric flow. For swept wings, two 
more bands are added both above and below the diagonal, the number of nonzero elements in a row 
is 13, and the system bandwidth is increased by 2*KXB. There are additional terms due to the 
wake that lie outside the banded system described above, and these are the same for both 
rectangular and swept planforms. The right-hand-side matrix, B, includes the wing boundary 
conditions. 

The general solution procedure developed for this report stores the coefficients in blocks using 
out-of-core storage. The basic banded set of coefficients is stored in A. The additional coefficients in 
A due to the wake are stored as part of B. The resulting system is solved using LU decomposition, 
and a Sherman-Morrison update procedure is then used to include the additional wake terms and 
thus obtain the complete solution (see ref. 24). 
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